#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <R.h>

int DSTP(double *zRS, double *zTS, 
          double *aURS, double *aUTS,
          double *aLRS, double *aLTS,
          double *ter,
          double *vT, double *vA, double *vSS, double *vRS2,
          double *svT, double *svA, double *svSS, double *svRS2,
          double *szRS, double *szTS, 
          double *ster,
          double *sRS, double *sTS,
          double *congruency,
          double *h, double *resp, double *rt,
          double *n, double *maxiter,
          int *rangeLow, int *rangeHigh, double *randomTable)
{
  double rhsRS,rhsTS,xRS,xTS,randNum,sampleVT,sampleVA,sampleVSS,sampleVRS2,currVRS,sampleTer,currV,time;
  int N,i,j,iter,Maxiter,blah,TS;
  
  N=(int) *n;
  Maxiter =(int) *maxiter;
  GetRNGstate();
  rhsRS=sqrt(*h)*(*sRS);
  rhsTS=sqrt(*h)*(*sTS);
  
  
  
  for (i=0;i<N;i++) {
    
    if (*szRS < 0.00001) {
      xRS=*zRS;
    } else {
      xRS = (*zRS) + ((rand()/(1.0 + RAND_MAX))*(*szRS)) - ((*szRS)/2);
    }
    
    if (*szTS < 0.00001) {
      xTS=*zTS;
    } else {
      xTS = (*zTS) + ((rand()/(1.0 + RAND_MAX))*(*szTS)) - ((*szTS)/2);
    }
    
    
    if (*svT < 0.00001) {
      sampleVT=*vT;
    } else {
      randNum = rand()/(1.0 + RAND_MAX);
      blah = (*rangeHigh) - (*rangeLow) + 1;
      blah = (randNum * blah) + (*rangeLow);
      randNum = randomTable[blah];
      
      sampleVT = (*vT) + ((*svT)*randNum);
    }
    
    if (*svA < 0.00001) {
      sampleVA=*vA;
    } else {
      randNum = rand()/(1.0 + RAND_MAX);
      blah = (*rangeHigh) - (*rangeLow) + 1;
      blah = (randNum * blah) + (*rangeLow);
      randNum = randomTable[blah];
      
      sampleVA = (*vA) + ((*svA)*randNum);
    }
    
    if (*svSS < 0.00001) {
      sampleVSS=*vSS;
    } else {
      randNum = rand()/(1.0 + RAND_MAX);
      blah = (*rangeHigh) - (*rangeLow) + 1;
      blah = (randNum * blah) + (*rangeLow);
      randNum = randomTable[blah];
      
      sampleVSS = (*vSS) + ((*svSS)*randNum);
    }
    
    if (*svRS2 < 0.00001) {
      sampleVRS2=*vRS2;
    } else {
      randNum = rand()/(1.0 + RAND_MAX);
      blah = (*rangeHigh) - (*rangeLow) + 1;
      blah = (randNum * blah) + (*rangeLow);
      randNum = randomTable[blah];
      
      sampleVRS2 = (*vRS2) + ((*svRS2)*randNum);
    }
    
    
    if (*ster < 0.00001) {
      sampleTer=*ter;
    } else {
      sampleTer = (*ter) + ((rand()/(1.0 + RAND_MAX))*(*ster)) - ((*ster)/2);
    }
    
    
    iter=0;
    resp[i]=(double) -1.0 ;
    TS=0;
    currV = sampleVT + (sampleVA*(*congruency));
    do 
    {
      iter = iter+1;
      
      if (TS == 0) {
        randNum = rand()/(1.0 + RAND_MAX);
        blah = (*rangeHigh) - (*rangeLow) + 1;
        blah = (randNum * blah) + (*rangeLow);
        randNum = randomTable[blah];
        
        xTS = xTS + (sampleVSS*(*h)) + (rhsTS*randNum);
        
        if (xTS>=*aUTS) {
          TS=(int) 1;
          currV = sampleVRS2;
        }
        if (xTS<=*aLTS) {
          TS=(int) 1;
          currV = sampleVRS2*(*congruency);
        }
      }
      
      
      randNum = rand()/(1.0 + RAND_MAX);
      blah = (*rangeHigh) - (*rangeLow) + 1;
      blah = (randNum * blah) + (*rangeLow);
      randNum = randomTable[blah];
      
      xRS = xRS + (currV*(*h)) + (rhsRS*randNum);
      
      if (xRS>=*aURS) {
        resp[i]=(double) 1.0 ; 
        break ;
      }
      if (xRS<=*aLRS) {
        resp[i]=(double) 2.0 ; 
        break ;
      }
    } while (iter<Maxiter) ; 
    rt[i]=(((double) iter)*(*h) - (*h)/((double) 2.0)) + sampleTer;
  }
  PutRNGstate();
}

